######################################################
# CBQ functionality
######################################################

library(rstan)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())


cbq <- function(N, # number of obs
                n_covariate, # number of covariates
                X, # covariates
                Y, # d.v
                qtl, # quantile to be estimated 1...9
                nchain = 2, # number of parallel chains
                niter = 1000, # number of interations
                thin = 1, # thinning
                seed = 1234567, # random seeds
                offset = 1e-20, # offset to enhance numerical stability
                max_treedepth = 10
                ){
  datlist = list(N=N,
                 Y = Y,
                 D = n_covariate,
                 X = X)
  pars = c("alpha","beta")
  ms_cbq_model = stan_model(paste("simulation/stan_new_models/cbq_binary_q",qtl,".stan",sep=""))
  ms_cbq_stan = sampling(ms_cbq_model, 
                         data=datlist,
                         pars=pars, 
                         seed=seed,
                         chains=nchain,
                         iter=niter,
                         thin=thin,
                         control = list(max_treedepth = max_treedepth))
  return(ms_cbq_stan)
}


pald <- function(x,p){
  out <- ifelse(x<0,p*exp(x*(1-p)), 1-(1-p)*exp(-x*(p)))
  return(out)
}